Thanks Delphine,
I have tested with float wmo=6902865, splitting the time series at cycle 80.

Here are the figures that have changes between intitial execution (left) and 4.5 run (right)

Figure 1

Figure 3

Figure 5

Figure 6

Suggested changes:

1. For compatibility with the previous version, graph_format = 'eps' should be kept as the default value
   l 37: graph_format = 'eps' ; % 'png' or 'eps'

2. l 302 & l 322 : change the position of the legend, add a box

legend([p1 p2 p3],{'2 x cal error','1 x cal error','1-1 profile fit'}, 'Location', 'Southoutside','Orientation','horizontal');
box on

3. l 407 : change ceil(max(TEMP)) by ceil(max(TEMP,[ ],'all'))

fl.yaxes = [floor(min(tlevels,[],'all')) ceil(max(tlevels,[],'all')) ceil(max(TEMP,[],'all'))];

Thanks Cécile for your comments. I've updated the Pull Request accordingly. Let me know if everything is fine now.

For Figure 5, if the upper most stable theta level is, say, around 12°C, and the lower most stable theta level is around 4°C (e.g. in NE Atlantic), then I would expect that most of the figure to be squeezed into the lower panel.
So it could be (~l 407 ):
fl.yaxes = [floor(min(min(TEMP,[],'all'),2)) ceil(min(max(tlevels,[],'all'),5)) ceil(min(max(TEMP,[],'all'),20))];

This ensure that the bottom panel is always a zoom on the deepest levels.
What do you think?

The change should also be made for figure 7 (~l.664)

I think this is a good idea; I've updated the code as suggested.